Mechanical analogue for cities

Motivated from the increasing need to develop a science-based, predictive understanding of the dynamics and response of cities when subjected to natural hazards, in this paper, we apply concepts from statistical mechanics and microrheology to develop mechanical analogues for cities with predictive capabilities. We envision a city to be a matrix where cell-phone users are driven by the city’s economy and other incentives while using the collection of its infrastructure networks in a similar way that thermally driven Brownian particles are moving within a complex viscoelastic material. Mean-square displacements of thousands of cell-phone users are computed from GPS location data to establish the creep compliance and the resulting impulse response function of a city. The derivation of these time-response functions allows the synthesis of simple mechanical analogues that model satisfactorily the city’s behaviour under normal conditions. Our study concentrates on predicting the response of cities to acute shocks (natural hazards) that are approximated with a rectangular pulse; and we show that the derived solid-like mechanical networks predict that cities revert immediately to their pre-event response suggesting an inherent resilience. Our findings are in remarkable good agreement with the recorded response of the Dallas metroplex following the February 2021 North American winter storm.

complex networks have been created and upgraded independently during different times with no apparent common design principles. Urbanization is an ever-growing, pressing challenge. As recently as 1950, only 30% of the world's population lived in cities. Today more than 50% is urbanized, and by 2050 more than 3/4 of the world's population is expected to live in cities [7,8]. As cities continue to grow, many of them along the coasts of continents which are prone to natural hazards, urban resilience has become a subject of interest in various disciplines [9][10][11][12][13]. The many challenges associated with rapid urban growth, community resilience and sustainable development are receiving increasing attention; yet, in practice, they have been treated as independent issues [14]. Linear per capita indicators that conflate distinct dynamic features pertinent to the dynamics specific to each city are typically used to assess these growing challenges; however, it appears these fail to yield direct measures of the impact that hazards can have on cities (such as earthquakes or hurricanes) in association with the need to quantify urban resilience and functional recovery [9,15].
Several parallels have been suggested between cities and other complex systems-from biological organisms [16][17][18][19] to insect colonies [20] and ecosystems [21]. Cities, while they appear as web-like structures, are much more than biological organisms or anthills, as they involve perpetual exchanges among their citizens who in their majority are productive, creative and innovative individuals, striving for upward mobility and social interaction. Accordingly, while the functionality of cities is served with the superposition of a variety of large-scale, complex networks (transportation, water, electricity, telecommunication) with no apparent common design principles, it appears that cities display a difficult-to-interpret self-organizing behaviour that leads to deterministic patterns that can be possibly described with an emerging mechanism yet to be identified [4,8,22]. Past studies [23][24][25] have uncovered that, in contrast to the random trajectories predicted by the initially proposed random walk models [26], human trajectories show a high degree of temporal and spatial regularity.
Part of the success of natural sciences hinges upon reductionism where analysis frameworks and methodologies are developed for predicting the emergent macroscopic behaviour of a system by monitoring and analysing the behaviour and/or properties of its constituents. For instance, in the emerging field of microrheology, the bulk frequency-and time-response functions of complex viscoelastic materials are inferred by monitoring the thermally driven Brownian motion of probe microspheres immersed within the viscoelastic material and subjected to the perpetual random forces from the collision of the molecules of the viscoelastic material [27][28][29][30][31][32][33][34][35][36][37].
Building on past studies on human mobility patterns [23][24][25] which uncovered that citizens exhibit characteristic travel distances in association with a significant probability to return to a few highly frequent locations-a behaviour that is reminiscent to particles in a potential well; in this work, we envision a city to be a matrix where people (cell-phone users) are driven by the city's economy and other personal and societal driving incentives while using the collection of its infrastructure networks in a similar way that thermally driven Brownian probe particles are moving within a complex viscoelastic material [27][28][29][30][31][32][33][34][35][36][37]. Mobile phones carried by the citizens of an urban centre offer rich data on human mobility patterns and we show that the ensemble averages of human GPS location data reveal information that characterizes the dynamics of the entire urban centre. Upon computing the mean-square displacements (MSD-ensemble averages) of the recorded time histories of the paths of a large number of citizens of a given city, we propose deterministic mechanical networks with a creep compliance that is proportional to the measured MSD. The response analysis of the proposed solid-like networks offers new insights on the mechanics of urban centres and their emergent inherent resilience within the context of engineering resilience [38].
Meerow et al. [13] identify several variations of the definition of urban resilience in the engineering, business/finance and social science literature. Our study concentrates on a single-state equilibrium (or a steady-state regime at equilibrium) which refers to the capacity of a system to revert to its postdisturbance equilibrium state [38]; therefore, adopts the generic definition of engineering resilience-that is the capacity of a system to recover its initial state and resume normal activity after a shock [39][40][41] and hinges upon the premise that a dependable indicator of the engineering resilience of a city [38] is whether the average mobility pattern of its citizens following an acute shock matches the average mobility pattern before the shock. Our work contributes to a quantitative, science-based, predictive understanding of the dynamics and response of cities to acute shocks [4,8] and concludes that large cities of average population density following an acute shock revert immediately to their initial steady-state, pre-event behaviour.

Ensemble averages of GPS locations of cell-phone users (IDs)
While the trace of the locations of each individual cell-phone ID shown in figures 1 and 2 is unique and every individual has a distinct purpose for generating the recorded trace, the collection of all traces of the royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943 anonymous cell-phone users in a given city is treated in this work as a random (stochastic) process. For each of the metroplexes examined in this work we processed longitude and latitude data (in degrees from the Greenwich meridian and from the equator) at various times for tens of thousands of anonymous cellphone users (IDs). The data were processed and the origin of every ID (home) was identified by extracting the most frequent occupied location by the ID during the time interval 01.00 to 04.00. With reference to figure 3, the origin = 'home' of every ID is established; the east-west (E-W, longitude) and north-south (N-S, latitude) displacement of each ID j at time t k is calculated from [42] and where R E = radius of Earth ≈ 6371 km and the displacement of ID j at time t k is r j ðt k Þ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi figure 3. For higher precision, one can use the haversine formula; yet the differences in the results from equations (2.1) and (2.2) are negligible given the relative small bounding areas of interest shown in figures 1 and 2. From equations (2.1) and (2.2), when ID j moves to the east of its origin (home), H j , its displacement projection, x j , is negative; whereas, when it moves to the west of its home, H j , its displacement projection, x j is positive. Similarly, when ID j moves to the north of its home, H j , its displacement projection, y j , is positive; whereas, when it moves to the south of its home, H j , its displacement projection, y j , is negative. Accordingly, for every city, the GPS location data of M IDs are organized as shown in table 1. For every time t k appearing in the first column of table 1 (say top of the hour), the location x j (t k ), y j (t k ) is the GPS location of ID j recorded closer to that time t k . Figure 4 plots time-histories of the E-W movement, x j (t), j ∈ {1, 2, … , M}, of selected IDs from the Dallas metroplex. Some IDs move frequently both to the east and to the west, other IDs move systematically to the east, other IDs move systematically to the west due to their daily routine (say travelling to their workplace). Accordingly, the time-averages from individual time-histories may differ drastically, therefore, the stochastic (random) process shown in figure 4 is clearly non-ergodic [43,44]. With reference to figure 4 for every time t k , we compute ensemble averages of the E-W movement 〈x(t k )〉 and the N-S movement 〈y(t k )〉 (mean values) for all the M available IDs . Displacement r j (t k ) of ID j at time t k from its 'home' H j . In statistical mechanics and microrheology, the symbol Δr j (t k ) rather than r j (t k ) is often used to express the distance of a probe from its 'origin'. In this work, the symbol Δ is dropped since the displacement r j is always measured from the fixed origin H j .  Similarly, for every time t k and a given time shift, τ, we compute the location ( position) autocorrelation functions of the E-W movement 〈x(t k )x(t k + τ)〉 and the N-S movement 〈y(t k )y(t k + τ)〉 for all the M available IDs in a city (ensemble).
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943 and y j ðt k Þy j ðt k þ tÞ: ð2:4bÞ   royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943 and 24 h. Clearly, for τ = 24 h, the location autocorrelation functions shown in figure 6 show a stronger correlation than when τ = 12 h, since after 24 h the cell-phone users are likely to be at the same location that they were at the day before. Accordingly, while the mean values of the locations (positions) shown in figure 5 fluctuate within only a couple of hundred metres, the strong fluctuations over time of the position autocorrelation functions shown in figure 6 indicate that the stochastic (random) process schematically illustrated in figure 4 does not even satisfy the notion of a weakly stationary process [43,44]. Accordingly, in this study, we work with ensemble averages (not time averages).

Theory
The analysis framework that we develop herein to construct a mechanical analogue for cities is inspired from major advances in microrheology in the mid-to late 1990s [27][28][29][30][31][32][33][34][35][36][37] which make possible the extraction of the bulk mechanical properties of complex viscoelastic materials by monitoring the thermally driven Brownian motion of probe microspheres suspended within the viscoelastic material. The phenomenon of Brownian motion was first explained in Einstein's celebrated 1905 paper [45] which examined the long-term response of Brownian microspheres with mass m and radius R suspended in a memoryless Newtonian fluid with viscosity η. Einstein's theory of Brownian motion predicts the long-term expression of the MSD, 〈r 2 (t)〉 = MSD, of the randomly moving microspheres (diffusive regime) where M is the number of suspended microspheres, while, r j (t) is the distance of microsphere j at time t from some origin. In equation (3.1), N ∈ {1, 2, 3} is the number of spatial dimensions, K B is Boltzmann's constant, T is the equilibrium temperature of the Newtonian fluid with viscosity η within which the M Brownian microspheres are immersed and D = K B T/6πRη is the time-independent diffusion coefficient. Equation (3.1) derived by Einstein in 1905, in association with Stokes' result that the drag force on a slowly moving sphere with velocity v, is 6πRηv [46], is of central interest in statistical mechanics since it relates the thermally driven ensemble average 〈r 2 (t)〉 from a stochastic process to the deterministic emergent property of the Newtonian fluid-that is its viscosity, η. The time-derivative of equation (3.1) gives The reader recognizes that the right-hand side term of equation (3.2) is a constant =NK B T/3πRη = 2ND; whereas, the left-hand side ð2=MÞ P M j¼1 j r j ðtÞðdr j ðtÞ=dtÞ j is zero when t = 0, since r j ð0Þ ¼ 0, 8j [ f1, 2, . . . , Mg. This inconsistency emerges because equation (3.1) is valid only in the longterm (diffusive regime at large times). At short timescales, when t < m/6πRη = τ, the Brownian motion of suspended particles is influenced by the inertia of the particle and the surrounding fluid (ballistic regime) and Einstein's 1905 'long-term' result offered by equation (3.1) was extended for all timescales by Uhlenbeck & Ornstein [47].

Generalization for Brownian motion within non-Newtonian fluids
In their seminal paper, Mason & Weitz [27] employed dynamic light scattering to measure the MSD 〈r 2 (t)〉 of probe Brownian particles immersed in a linear viscoelastic fluid and related it to the complex dynamic modulus, G ve ðsÞ, of the viscoelastic fluid within which the probe particles are immersed. Given the viscoelastic behaviour of the complex fluid, the motion of a probe Brownian particle is described by the generalized Langevin equation [48,49] m dvðtÞ dt ¼ À where m is the mass of the Brownian particle, v(t) is the particle velocity and f R (t) are the random forces royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943 acting on the particle from the collisions of the fluid molecules on the Brownian particle. The integral in equation (3.4) represents the drag force on the particle as it moves through the viscoelastic fluid and accounts for the fading memory of this drag due to the elasticity of the viscoelastic fluid. Upon transforming equation (3.4) in the Laplace domain, Mason & Weitz [27] reached the following result for the Laplace transform of the MSD: where s is the Laplace variable. In deriving equation (3.5) a key assumption was adopted-that the Stokes drag coefficient on a sphere moving slowly in a memoryless, Newtonian viscous fluid, ζ = 6πRη [46], can be generalized to relate the impedance of the Brownian particle, ZðsÞ ¼ LfzðtÞg¼ Ð 1 0 À zðtÞ e Àst dt to the complex dynamic viscosity of the viscoelastic material, h ve ðsÞ ¼ G ve ðsÞ=s. Accordingly, Mason & Weitz [27] made the physically motivated assumption that In a recent publication, Makris [50] showed that the time response function appearing on the right-hand side of equation (3.3) is the creep compliance (retardation function) of a linear network where a dashpot with viscosity η is connected in parallel with an inerter with distributed inertance m R = m/6πR. Building on the work of [47,51] in association with the work of Mason & Weitz [27], Makris [50] uncovered a viscous-viscoelastic correspondence principle for Brownian motion which states that the MSD 〈r 2 (t)〉 of Brownian particles (microspheres) with mass m and radius R suspended in some linear, isotropic viscoelastic material when subjected to the random forces from the collisions of the molecules of the viscoelastic material, is where J(t) is the creep compliance (strain history due to a unit-step stress) of a viscoelastic network that is a parallel connection of the linear viscoelastic material within which the Brownian particles are immersed and an inerter with distributed inertance m R = m/6πR. Equation (3.7) also holds for the case where the density of the fluid surrounding the Brownian microparticles is appreciable [52] and in this case, in addition to the Stokes viscous drag, the Brownian motion of the immersed microparticles develops the Boussinesq hydrodynamic memory [53]. Equation (3.7), which is the equivalent of equation (3.5) in the time domain, generalizes the main result of equation (3.1) and relates the thermally driven ensemble average 〈r 2 (t)〉 = MSD from a stochastic (random) process to the deterministic creep compliance (retardation function), J(t) of a mechanical network. In this work, we build upon equation (3.7) in an effort to propose a corresponding relation for cities, where 〈r 2 (t)〉 is the MSD of cell-phone users as computed by averaging over a city ensemble of GPS location traces, J(t) is the creep compliance (displacement history due to a unit-step force) of a mechanical model to be identified and W is a proportionality constant having the units of energy (force × length). The motion of people in a city as recorded with their GPS locations happens in a matterless environment (the environment of a city which does not exhibit mechanical elasticity, dissipation or inertia); while it is driven by conceptual incentives; therefore, we do not have a Langevin equation to describe the motion of people in cities as we have equation (3.4) that describes the random motion of Brownian particles immersed in some linear, isotropic viscoelastic material. Furthermore, the statistical analysis of the GPS location data processed in this study yields that ensemble averages as illustrated in figure 4 do not even support the notion of a weak stationary random process [43,44]. Nevertheless, in view of the overall phenomenological similarities between the motion of people in a city (figures 1 and 2) and the motion of probe Brownian microparticles in a viscoelastic material in association with the remarkable success of equations (3.5) or (3.7) in microrheology, we adopt equation (3.8) in an effort to develop a mechanical model for cities with potential engineering significance. Figure 7 plots royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943 the time history of the MSD = 〈r 2 (t)〉  Figure 7 reveals that the MSD for the Dallas metroplex exhibits a remarkable periodicity with people moving more during the weekdays and less on Sundays which are marked with green dots. At the same time, less people are staying at home during the weekends than during the weekdays. The MSD from the Dallas metroplex exhibits a notable suppression during the major North American winter storm that happened during the third week of February 2021 and ended sometime towards the end of that week before Sunday, 21 February 2021. Over 4 million people lost power due to the storm which resulted to a state-wide power crisis in Texas causing major damage to the Texas power grid [54,55]. Regardless of its unprecedented intensity and resulting disruption from the winter storm, figure 7 reveals that immediately after the cold winter storm, the city of Dallas resumes normal activity with the MSD of its citizens reverting to its normal oscillation pattern as recorded before the storm, suggesting that the Dallas metroplex exhibits a great degree of engineering resilience [38].
The MSD from the San Francisco (SF) Bay Area for the same period (1 February-31 March 2021) is less ordered that the MSD from the Dallas metroplex; nevertheless the same pattern of less people staying at home during the weekends than during the weekdays is observed. Furthermore, the average MSD from the SF Bay Area covers less area (on an average 60 km 2 ) when compared with  The smaller values of the average MSD of the SF Bay Area when compared with that from the Dallas metroplex is partly due to the presence of the water boundaries (Pacific Ocean and SF Bay) in association with the presence of the Santa Cruz mountains to the southwest (figure 2). At the same time, the SF Bay Area is served with pronounced traffic corridors which are connected with the major bridges over the Bay Area leading to different traffic patterns than these of the landlocked Dallas metroplex. Nevertheless, the recorded MSD for the San Francisco Bay Area shown in figure 8 and its satisfactorily match with the proposed amplitude-modulation model shown in figure 13 suggest that our mechanical model may be also promising for cities with pronounced traffic corridors and traffic patterns other than these of the Dallas metroplex. The cell-phone mobility data that we have access to at this time go back to 2019 only; and for this 4-year period, there has not been any natural-hazard (earthquake) that struck the San Francisco Bay Area to confirm the predictions of the proposed mechanical model for such cities. Figures 7 and 8 show that the MSD, 〈r 2 (t)〉, of cell-phone users as computed by averaging over an ensemble of citizens from major urban areas, is essentially an elevated cosine with a daily period where the amplitude of the cosine oscillation is modulated with a weekly period. Accordingly in view of equation (3.8), in this section we are in search of mechanical networks that their creep compliance, J(t) (displacement history due to a unit-step force) can approximate the time history of the MSD, 〈r 2 (t)〉 as computed from ensemble averages and shown in figures 7 and 8.

Candidate mechanical analogues for cities
Given that the MSD, 〈r 2 (t)〉, shown in figures 7 and 8 are bounded functions, in association with equation (3.8) we are in search of a solid-like network so that its creep compliance, J(t), is a bounded time response function. Furthermore, the perpetual undamped oscillations of the computed 〈r 2 (t)〉 with no visible rise-time imply that the candidate mechanical network shall be undamped. Finally, given that the cosine function shown in figures 7 and 8 is elevated above the horizontal axis, the candidate network needs to yield a creep compliance function that is more elaborate than the elementary inertoelastic solid (the spring-inerter parallel connection) [56,57].

The three-parameter inertoelastic solid
The simplest mechanical network that meets the above requirements is a parallel connection of an elastic spring with an inertoelastic fluid (spring-inerter in-series connection) as shown in figure 9. The total force F(t) = F 1 (t) + F 2 (t) from the linear network shown in figure 9 is the summation of the force output from the linear spring with elastic constant k 1 ,  FðtÞ HðsÞ ¼ uðsÞ The Laplace transform of the creep compliance, J(t) is the complex creep function CðsÞ ¼ LfJðtÞg¼ Ð 1 0 JðtÞ e Àst dt ¼ HðsÞ=s, [58][59][60][61]. Upon dividing equation (4.5) with the Laplace variable s, we obtain CðsÞ ¼ HðsÞ The inverse Laplace transform of equation (4.6) offers the creep compliance of the three-parameter inertoelastic solid shown in figure 9.
where U(t − 0) is the Heaviside unit-step function at the time origin [62] and parameter b ¼ v 2 R2 =v 2 R1 ¼ k 2 =k 1 . The expression of the creep compliance, J(t) offered by equation (4.7) is an elevated cosine function with frequency Substitution of equation (4.7) into equation (3.8), which is the central equation for this study, the elementary three-parameter inertoelastic solid mechanical analogue predicts that  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943 where the proportionality constant W/k 1 = A 0 has units of area (km 2 ). Accordingly, given that is the daily frequency (ω R = π/12 rad h −1 ), only two parameters of the three-parameter inertoelastic solid need to be calibrated from the recorded GPS location data;

An amplitude-modulation model
The time history of the MSD 〈r 2 (t)〉 of the Dallas metroplex IDs shown in figure 7 exhibits a distinguishable weekly modulation (message signal) in addition to the daily oscillations (carrier wave), C(t) = Ccos (ω R t) with ω R = 2π/24 h = π/12 rad h −1 . By expressing the weekly modulating message signal [63,64] an improved mathematical expression to equation (4.8) that better approximates the Dallas metroplex MSD shown in figure 7 is where C and γ are dimensionless model parameters. Upon using the trigonometric identity cosðv R tÞ sin (ðv R =7Þt) ¼ 1 2 ½sinð 8 7 v R tÞ À sinð 6 7 v R tÞ, equation (4.11) assumes the form

Response of the candidate models to a rectangular pulse with duration T e
We return now to the mechanical model depicted in figure 9 which is the simplest mechanical analogue for a city since its creep compliance J(t) given by equation (4.7) offers according to equation (3.8) a MSD 〈r 2 (t)〉 that matches to a satisfactory extent the recorded MSD as results from the GPS location data of several thousand users (figures 10 and 11). The main motivation of this study is to develop quantitative models for cities that have some predictive capability. Our response analysis starts with an elementary rectangular force pulse of amplitude F o and duration T e which is expressed as ð5:1Þ where U(t − ξ) is the unit-amplitude Heaviside step function that initiates at time t = ξ [62]. The rectangular pulse with duration T e expressed by equation (5.1) is suitable to express shocks from natural hazards (earthquakes, hurricanes and cold storms) which start and end abruptly and are system-wide [65]. By contrast, Gaussian functions with long tails are less suitable to approximate natural hazards which impact cities suddenly with little or zero loading build-up time. Pandemics and economic crises are stressing cities during larger timescales during which the city may develop mechanisms to readjust [13,66]. Therefore, for such prolonged crisis, different loading patterns than the rectangular pulse described with equation (5.1) may be more suitable.
The response of the three-parameter inertoelastic model shown in figure 9, which is described with the constitutive equation (4.3), to the forcing function given by equation (5.1) can be computed with equation (4.5), which gives uðsÞ ¼ HðsÞFðsÞ, in which F(s) is the Laplace transform of F(t) given by equation (5.1) and HðsÞ is the complex dynamic flexibility, given by equation (4.5). Accordingly, under a rectangular pulse-load s HðsÞ e ÀsTe : ð5:3Þ

Response of the three-parameter inertoelastic solid to a rectangular pulse
Recognizing that HðsÞ=s ¼ CðsÞ is the complex creep function expressed by equation (4.6) and after expanding the first term in the bracket of the right-hand side of equation (4.5) in partial fractions the inverse Laplace transform of equation (5.3) gives cos½v R ðt À T e Þ : ð5:5Þ Upon replacing in equation (5.5) the expression of J(t) given by equation (4.7) and after performing the remaining inverse Laplace transforms appearing in equation (5.5), the response of the three-parameter inertoelastic solid shown in figure 9 when subjected to the rectangular pulse loading given by equation (5.1) is Clearly, when the duration of the acute shock (rectangular pulse) is an integer multiple of a day (T e = n(24 h) = n(24 h/2π)2π = n(2π/ω R )), the last two cosine terms in equation (5.6) cancel and the displacement response reduces to u(t) = ( On the other hand, when the duration of the rectangular pulse is not an integer multiple of a day (say T e = 138 hours), the displacement response of a city, u(t), upon being subjected to the rectangular pulse exhibits perpetual oscillations ( figure 14c). Consider now a city operating under normal conditions, therefore its MSD 〈r 2 (t)〉 of its citizens is expressed with equation (4.8). For the simplest admissible mechanical analogue shown in figure 9, the corresponding normalized MSD k 1 〈r 2 (t)〉/W = 〈r 2 (t)〉/A o given by equation (4.8) is plotted at the top of figure 14a. While the city operates under normal conditions, it is suddenly subjected to an acute shock at some time t 0 = 0 that is described with the rectangular pulse given by equation ( in figure 14c. Given that in this work, we consider the city to be a linear network, the overall response of the city subjected to the acute shock (rectangular pulse) is the superposition of its steady-state response under normal conditions given by equations (4.7) or (4.8) and shown in figure 14a (response to a unit step loading at the distant past that produces the creep compliance, J(t)) and its response due to the rectangular pulse at time t 0 given by equation (5.6) and shown in figure 14c. Accordingly, the 'postevent' response of a city as is described with the three-parameter inertoelastic solid is The minus sign (rather than a plus sign) in the left-hand side of equation (5.7) is because the acute shock is expected to suppress the normal activity of the city as expressed by equations (4.7) or (4.8).
3 city under normal conditions described with equations (4.7) or (4.8) acute shock expressed in this study with a rectangular pulse with duration T e : response of the mechanical analogue for a city due to acute shock alone equation ( royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943 The right-hand side of equation (5.7) that results after subtracting equation (5.6) from equation (4.7) is precisely of the same form as equation (4.7) and uncovers the remarkable result that upon the expiration of the acute shock (rectangular pulse), the mechanical model shown in figure 9 and expressed by equation (4.3) predicts that the city reverts immediately to its normal activity (same 〈r 2 (t)〉 as before the acute shock), implying that cities when subjected to a rectangular pulse of any duration, T e are inherently and invariably resilient. Figure 14d plots the pre-and post-event response as predicted by equation (5.7) and compares the analytical prediction with the recorded response of the Dallas metroplex following the 2021 North America cold storm during the third week of February 2021 for which we have used a pulse excitation duration T e = 138 h. The result that cities of average population density such as the Dallas metroplex examined in this study are inherently resilient to acute shocks as indicated by equation (5.7) is also illustrated in acute shock expressed in this study with a rectangular pulse with duration T e : response of the mechanical analogue for a city due to acute shock alone equation ( royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943 figure 15 for the case where the duration of the rectangular shock is an integer multiple of a day (T e = 144 h) as shown in figure 15b. Given the linearity of the candidate constitutive model for a city expressed by equation (4.3) and illustrated in figure 9, an alternative way to reach equation (5.6) is through the convolution integral where F(ξ) is the time history of the acute shock-in this analysis the rectangular pulse given by equation (5.1)-and h(t) is the impulse response function [67][68][69][70][71] that is defined as the resulting displacement history u(t), due to an impulsive force input F(t) = δ(t − ξ) with ξ < t. Given that δ(t − ξ) = dU(t − ξ)/dt, the impulse response function of a linear system h(t) = dJ(t)/dt where J(t) is its creep compliance (displacement history due to a unit step force). Accordingly, the impulse response function of the three-parameter inertoelastic solid described with equation (4.3) and illustrated by figure 7 is where J(t) is given by equation (4.7) and δ(t − 0) is the Dirac delta function [62]. Substitution of the expression of the impulse response function given by equation (5.9) into equation (5.8) after using the forcing function F(ξ) given by equation (5.1) yields Using the symmetry of the Dirac delta function δ(t − ξ) = δ(ξ − t) [62], the first integral appearing within the brackets of equation (5.10) gives and after expanding the sine function in the second integral, equation (5.10) assumes the form Upon evaluating the cosine and sine integrals in the right-hand side of equation (5.12), the displacement response history offered by equation (5.6) is recovered. This second approach, where the displacement response is evaluated with the convolution integral given by equation (5.8), is employed to calculate the response from the amplitude-modulation model for which we do not have a complex frequency response function, HðsÞ at this time.

Response of the amplitude-modulation model to a rectangular pulse
For a city that the MSD of its citizens 〈r 2 (t)〉 is best described with the amplitude-modulation model expressed with equation (4.11) the corresponding creep compliance is Again, when the duration of the rectangular pulse is not an integer multiple of a day (say T e = 138 h), the displacement response of a city u(t) described with the amplitude modulation model exhibits perpetual oscillations as shown in figure 16c. Consider now a city operating under normal conditions and the MSD 〈r 2 (t)〉 of its citizens is best described by the amplitude-modulation model given by equation (4.11) and plotted at the top of figure 16a. While the city operates under normal conditions, at some time t 0 = 0 is suddenly subjected royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943 to an acute shock that is described with a rectangular pulse with duration T e as shown in figure 16b. In the response analysis shown in figure 16, we have selected a duration of the pulse-excitation to be a noninteger multiple of a day (T e = 138 h = 5 days + 18 h) so that the city response u(t) as expressed by equation (5.15) exhibits oscillations. The overall response upon the city is subjected to the acute shock (post-event response) is the superposition of its steady-state response under normal conditions given by equation (4.11) and shown in figure 16a and its response due to the rectangular pulse given by equation (5.15) and shown in figure 16c. Accordingly, the 'post-event' response of a city described with the amplitude-modulation model is acute shock expressed in this study with a rectangular pulse with duration T e : response of the city-model due to acute shock alone: equation royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943 The right-hand side of equation (5.16), that results after subtracting equation (5.15) from equation (4.11), is precisely of the same form as equation (4.11) and confirms our previous finding that, regardless of the sophistication of the linear mathematical model used to approximate the MSD of the citizens of a city, upon the expiration of the acute shock (rectangular pulse), the city reverts immediately to its normal pre-event activities (same 〈r 2 (t)〉 as before the acute shock) implying that cities when subjected to a rectangular pulse are inherently resilient within the context of engineering resilience [38].
The remarkable result-that cities of average population density such as the Dallas metroplex study are inherently resilient to acute shocks, regardless the sophistication of the linear solid-like model that approximates the recorded MSD-is also illustrated in figure 17 for the case where the duration of the rectangular pulse is an integer multiple of a day as shown in figure 17b. For the amplitudemodulation model, even when the duration of the rectangular shock is an integer number of days (T e = 144 h = 6 days) the sine functions appearing in equation (5.15) do not cancel because of the multiplication factors 8/7 and 6/7, resulting to the small undulations that override the rectangular pulse shown in figure 17c. The post-event predictions of our mathematical models are in remarkably good agreement with the recorded response of the Dallas metroplex following the 2021 North American winter storm during which over 4 million people lost power due to a resulting state-wide power crisis in Texas causing unforeseen damage to its power grid [54,55].

Accurate estimation of the recovery time of cities by synchronizing the predicted and recorded post-event responses
While our mechanical model predicts with remarkable accuracy the recorded post-event response, it is not capable to capture the local, transient response of a city during the duration of the natural hazard. Our main motivation for developing a mechanical model for cities is for estimating the recovery time of the city when struck by a natural hazard-that is how much time the city needs to revert to its initial steady-state behaviour and resume its normal pre-event activities. Our response analysis framework offers an accurate method for estimating the recovery time of cities when struck by natural hazards by merely taking the duration of the acute pulse, T e , long enough so that the predicted postevent response synchronizes with the recorded post-event response as shown in figure 17d (or figure  15d ). More specifically, figure 16d shows that when we used as duration of the rectangular pulse, T e = 138 h, the predicted MSD of the Dallas metroplex following the storm is identical to the pre-event MSD; yet the predicted post-event MSD is slightly out of phase (6 hours) with the recorded MSD. By contrast, figure 17d shows the predicted post-event response of the Dallas metroplex when the duration of excitation of the rectangular pulse is T e = 6 days = 144 h, synchronizes with the recorded post-event MSD, indicating that the Dallas metroplex recovered essentially immediately from the February 2021 North American winter storm. In the general case, the recovery time, T rec , is the difference between T e (duration of the rectangular pulse) and the actual duration of the natural hazard, T nh , as documented by the federal and state agencies-that is T rec = T e − T nh . As an example, for the Dallas metroplex, when struck by the February 2021 North American winter storm, figure 17d or figure 15d suggest that the recovery time is essentially zero (T rec = T e − T nh = 0).

Conclusion
Following the ever-growing number of qualitative studies on urban resilience [13] in association with current trends for developing a quantitative, science-based, predictive framework for the behaviour/response analysis of cities [4,8,25]; in this work, we employ concepts from statistical mechanics and microrheology to develop a mechanical analogue for cities. Our study concentrates on a single-state equilibrium as defined by [38] which is the capacity of a system to revert to its post-disturbance equilibrium state-also known as engineering resilience-and hinges upon the premise that a dependable indicator of the engineering resilience of a city is whether the average mobility pattern of its citizens following an acute shock matches the average mobility pattern before the shock. We built on the remarkable success of microrheology [27][28][29][30][31][32][33][34][35][36][37] in which the macroscopic frequency and time-response functions of complex viscoelastic materials are extracted by monitoring the thermally driven Brownian motion of probe microparticles suspended within the viscoelastic material, and we computed the MSD (ensemble averages) of thousands of cell-phone users from GPS location data to royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 220943 establish the creep compliance and the resulting impulse response function of a city. The derivation of these time-response functions allows for the synthesis of simple solid-like mechanical analogues with increasing sophistication which are shown to have engineering significance.
Arguments from physical modelling led to the synthesis of the simplest mechanical analogue for a city, which is the three-parameter inertoelastic solid that is a parallel connection of an elastic spring with an inertoelastic fluid (a spring-inerter in-series connection). The three-parameter inertoelastic solid yields a creep compliance which is essentially an elevated cosine that captures the daily oscillations of the MSD of the citizens of a city in association that a city 'never sleeps'.
Recorded GPS location data from the Dallas metroplex and the San Francisco Bay Area indicate that in addition to the daily oscillations, the MSD of cities exhibits a weekly amplitude-modulation with people moving more during the weekdays and less on Sundays; while, less people are staying at home during the weekends than during the weekdays. These observations motivated the synthesis and response analysis of a solid-like, amplitude-modulation model which captures satisfactorily both the daily and weekly periodicity. Given their linear structure, both the three-parameter inertoelastic model and the more sophisticated amplitude-modulation solid-like model predict that, when a city is subjected to an acute shock that was modelled with a rectangular pulse with arbitrary finite duration, it recovers immediately its initial state ( pre-event response), suggesting that cities of average population density, such as the Dallas metroplex examined herein, are inherently resilient within the context of engineering resilience [38]. The mathematical predictions from the mechanical analogues developed in this study confirm similar findings by other investigators, reached with different analysis methodologies, that cities are remarkably robust [8].
Data accessibility. Recorded GPS location datasets for all related graphs reported in this paper are publicly accessible in the Dryad Digital Repository: https://doi.org/doi:10.5061/dryad.1c59zw3zg [72].